C
C       $Date: 1996/03/29 18:30:00 $
C       $Revision: 2.3 $
C
      SUBROUTINE PCG2AL(ISUM,LENX,LCV,LCSS,LCP,LCCD,LCHCHG,LCLHCH,
     1       LCRCHG,LCLRCH,MXITER,ITER1,NCOL,NROW,NLAY,IN,IOUT,NPCOND,
     A       LCIT1)
C
C-----VERSION 0002 01MAY1989 PCG2AL
C
C     ******************************************************************
C     ALLOCATE STORAGE IN THE X ARRAY FOR PCG ARRAYS
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
C     ------------------------------------------------------------------
C
C-------PRINT A MESSAGE IDENTIFYING PCG PACKAGE
      WRITE(IOUT,1)
    1 FORMAT(1H0,'PCG2 -- CONJUGATE GRADIENT SOLUTION PACKAGE'
     1,', VERSION 2.1, 6/1/95')
C
C-------READ AND PRINT MXITER,ITER1 AND NPCOND
      READ(IN,2) MXITER,ITER1,NPCOND
    2 FORMAT(3I10)
      WRITE(IOUT,3) MXITER,ITER1,NPCOND
    3 FORMAT(' MAXIMUM OF',I4,' CALLS OF SOLUTION ROUTINE'/
     1       ,' MAXIMUM OF',I4,' INTERNAL ITERATIONS PER '
     2       ,'CALL TO SOLUTION ROUTINE'/
     3       ,' MATRIX PRECONDITIONING TYPE :',I5)
C
C-------ALLOCATE SPACE FOR THE PCG ARRAYS
      ISOLD=ISUM
      NRC=NROW*NCOL
      ISIZ=NRC*NLAY
      LCV=ISUM
      ISUM=ISUM+ISIZ
      LCSS=ISUM
      ISUM=ISUM+ISIZ
      LCP=ISUM
      ISUM=ISUM+ISIZ
      LCCD=ISUM
      IF(NPCOND.NE.2) ISUM=ISUM+ISIZ
      LCHCHG=ISUM
      ISUM=ISUM+MXITER*ITER1
      LCLHCH=ISUM
      ISUM=ISUM+3*MXITER*ITER1
      LCRCHG=ISUM
      ISUM=ISUM+MXITER*ITER1
      LCLRCH=ISUM
      ISUM=ISUM+3*MXITER*ITER1
      LCIT1=ISUM
      ISUM=ISUM+MXITER*ITER1
C
C-------CALCULATE AND PRINT THE SPACE USED IN THE X ARRAY
      ICG=ISUM-ISOLD
      WRITE(IOUT,4) ICG
    4 FORMAT(1X,I7,' ELEMENTS IN X ARRAY ARE USED BY PCG')
      ISUM1=ISUM-1
      WRITE(IOUT,5) ISUM1,LENX
    5 FORMAT(1X,I7,' ELEMENTS OF X ARRAY USED OUT OF',I7)
      IF(ISUM1.GT.LENX) WRITE(IOUT,6)
    6 FORMAT(1X,'   ***X ARRAY MUST BE DIMENSIONED LARGER***')
C
C-------RETURN
      RETURN
      END
      SUBROUTINE PCG2RP(MXITER,ITER1,HCLOSE,RCLOSE,NPCOND,NBPOL,
     1                RELAX,IPRPCG,IN,IOUT,MUTPCG,NITER,IT1,DAMP)
C
C-----VERSION 0002 01MAY1989 PCG2RP
C     01SEPT1990 IPCGCD OMITTED; NITER ADDED
C     01JUNE1995 DAMP ADDED
C
C     ******************************************************************
C     READ DATA FOR PCG
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      DIMENSION IT1(MXITER*ITER1)
C     ------------------------------------------------------------------
C
C-------READ HCLOSE,RCLOSE,RELAX,NBPOL,IPRPCG,MUTPCG
      READ(IN,1) HCLOSE,RCLOSE,RELAX,NBPOL,IPRPCG,MUTPCG,DAMP
    1 FORMAT(3F10.0,3I10,F10.0)
C
C-------PRINT MXITER,ITER1,NPCOND,HCLOSE,RCLOSE,RELAX,NBPOL,IPRPCG,
C-------MUTPCG,DAMP
      WRITE(IOUT,100)
  100 FORMAT(1H0,///57X,'SOLUTION BY THE CONJUGATE-GRADIENT METHOD'
     1/57X,43('-'))
      WRITE(IOUT,115) MXITER
  115 FORMAT(1H0,38X,'MAXIMUM NUMBER OF CALLS TO PCG ROUTINE =',I9)
      WRITE(IOUT,120) ITER1
  120 FORMAT(1H ,42X,'MAXIMUM ITERATIONS PER CALL TO PCG =',I9)
      WRITE(IOUT,122) NPCOND
  122 FORMAT(1H ,49X,'MATRIX PRECONDITIONING TYPE =',I9)
      IF(NPCOND.EQ.2) WRITE(IOUT,123)
  123 FORMAT(1H ,58X,'THE MATRIX WILL BE SCALED')
      WRITE(IOUT,124) RELAX,NBPOL
  124 FORMAT(1H ,26X,'RELAXATION FACTOR (ONLY USED WITH',
     1' PRECOND. TYPE 1) =',E15.5,/,
     2 1H ,19X,'PARAMETER OF POLYMOMIAL PRECOND.'
     3 ,' = 2 (2) OR IS CALCULATED :',I9)
      WRITE(IOUT,125) HCLOSE
  125 FORMAT(1H ,43X,'HEAD CHANGE CRITERION FOR CLOSURE =',E15.5)
      WRITE(IOUT,127) RCLOSE
  127 FORMAT(1H ,39X,'RESIDUAL CHANGE CRITERION FOR CLOSURE =',E15.5)
      IF(IPRPCG.LE.0)IPRPCG=999
      WRITE(IOUT,130) IPRPCG,MUTPCG
  130 FORMAT(1H ,30X,'PCG HEAD AND RESIDUAL CHANGE PRINTOUT INTERVAL ='
     1,I9,/,1H ,23X,'PRINTING FROM SOLVER IS LIMITED(1) OR ',
     2'SUPPRESSED (>1) =',I9)
      IF(DAMP.LE.0.0) DAMP=1.0
      WRITE(IOUT,135) DAMP
 135  FORMAT(1H ,59X,'DAMPING PARAMETER =',E15.5)
      NITER=0
C
      RETURN
      END
      SUBROUTINE PCG2AP(HNEW,IBOUND,CR,CC,CV,HCOF,RHS,V,SS,P,CD,
     1      HCHG,LHCH,RCHG,LRCH,KITER,NITER,HCLOSE,RCLOSE,ICNVG,
     2      KSTP,KPER,IPRPCG,MXITER,ITER1,NPCOND,NBPOL,NSTP,NCOL,NROW,
     3      NLAY,NODES,RELAX,IOUT,MUTPCG,IU,IP,SN,SP,SR,IT1,DAMP)
C-----VERSION 0002 01MAY1989 PCG2AP
C     01JULY1990 COMMENT STATEMENTS ADDED AND MODIFIED
C     01SEPT1990 IPCGCD OMITTED; STATEMENT 590 ADDED
C     27SEPT1990 STATEMENT IN DO 155 LOOP CHANGED
C     01SEPT1991 ADDED STATEMENTS RELATED TO SENSITIVITY CALCULATIONS
C                AT THE END OF THE 115 LOOP.  CHANGED THE 510 FORMAT
C                STATEMENT AND THE PRECEDING IF STATEMENT
C     20MAR1992  CHANGED 510 FORMAT STATEMENT; OMITTED 2 LINES IN DO 160
C                LOOP
C     01MAY1993  ADDED DEL TO CALCULATION OF THE CHOLESKY DIAGONAL
C     15JUNE1993 MADE CELLS SURROUNDED BY DRY CELLS INACTIVE
C     01JUNE1995 ADDED DAMP
C
C     ******************************************************************
C     SOLUTION BY THE CONJUGATE GRADIENT METHOD -
C                                          UP TO ITER1 ITERATIONS
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      PARAMETER (DZERO=0.D0,DONE=1.D0)
      DOUBLE PRECISION HNEW,HHCOF,RRHS,DEL
      DOUBLE PRECISION Z,B,D,E,F,H,S,ALPHA
      DOUBLE PRECISION ZHNEW,BHNEW,DHNEW,FHNEW,HHNEW,SHNEW
      DOUBLE PRECISION SRNEW,SROLD,SSCR,SSCC,SSCV,VCC,VCR,VCV
      DOUBLE PRECISION CDCC,CDCR,CDCV
      DOUBLE PRECISION PN,VN,HCHGN,RCHGN,PAP
      DOUBLE PRECISION FCC,FCR,FCV,FV
C
      DIMENSION HNEW(NODES), IBOUND(NODES), CR(NODES), CC(NODES),
     1  CV(NODES), HCOF(NODES), RHS(NODES),IT1(MXITER*ITER1),
     2  V(NODES), SS(NODES), P(NODES), CD(NODES), HCHG(MXITER*ITER1),
     3  LHCH(3,MXITER*ITER1), RCHG(MXITER*ITER1), LRCH(3,MXITER*ITER1)
C     ------------------------------------------------------------------
C
C      IF(NITER.EQ.0) THEN
C        WRITE(IOUT,895)
C        WRITE(IOUT,900) (I,CC(I),CR(I),CV(I),HCOF(I),
C     1     RHS(I),HNEW(I),IBOUND(I),I=1,50)
C      ENDIF
  895 FORMAT('    I',5X,'CC',10X,'CR',10X,'CV',9X,'HCOF',8X,'RHS',9X,
     1   'HNEW',4X,'IBOUND (IBND=0 NOT PRINTED)')
  900 FORMAT(I5,5E12.3,D12.3,I5)
C-------ASSIGN VARIABLE EQUAL TO THE NUMBER OF CELLS IN ONE LAYER
      NRC=NROW*NCOL
C-------INITIALIZE VARIABLES USED TO CALCULATE ITERATION PARAMETERS
      SRNEW=DZERO
      BPOLY=0.
      IF(NPCOND.NE.1) RELAX=1.
      NORM=0
      IF(NPCOND.EQ.2) NORM=1
C-------INITIALIZE VARIABLE USED TO TEST FOR NEGATIVE CHOLESKY DIAGONAL
      CD1=0.
C------CLEAR PCG WORK ARRAYS.
      DO 100 N=1,NODES
      SS(N)=0.
      P(N)=0.
  100 V(N)=0.
C------FOR NPCOND=1, INITIALIZE CHOLESKY DIAGONAL
      IF(NPCOND.EQ.1) THEN
        DO 105 N=1,NODES
  105   CD(N)=0.
      ENDIF
C
C------CALCULATE THE RESIDUAL. IF NORM=1, CALCULATE THE DIAGONALS OF
C------THE A MATRIX,AND STORE THEM IN HCOF.
      DO 115 K=1,NLAY
      DO 115 I=1,NROW
      DO 115 J=1,NCOL
C
C-------CALCULATE 1 DIMENSIONAL SUBSCRIPT OF CURRENT CELL AND
C-------SKIP CALCULATIONS IF CELL IS INACTIVE
      N=J+(I-1)*NCOL+(K-1)*NRC
      IF(IBOUND(N).EQ.0) THEN
        CC(N)=0.
        CR(N)=0.
        IF(N.LE.(NODES-NRC)) CV(N)=0.
        IF(N.GE.2) CR(N-1)=0.
        IF(N.GE.NCOL+1) CC(N-NCOL)=0.
        IF(N.LE.(NODES-NRC).AND.N.GE.NRC+1) CV(N-NRC)=0.
        HCOF(N)=0.
        RHS(N)=0.
        GO TO 115
      ENDIF
C
C-------CALCULATE 1 DIMENSIONAL SUBSCRIPTS FOR LOCATING THE 6
C-------SURROUNDING CELLS
      NRN=N+NCOL
      NRL=N-NCOL
      NCN=N+1
      NCL=N-1
      NLN=N+NRC
      NLL=N-NRC
C
C-------CALCULATE 1 DIMENSIONAL SUBSCRIPTS FOR CONDUCTANCE TO THE 6
C-------SURROUNDING CELLS.
      NCF=N
      NCD=N-1
      NRB=N-NCOL
      NRH=N
      NLS=N
      NLZ=N-NRC
C
C-----GET CONDUCTANCES TO NEIGHBORING CELLS
C-------NEIGHBOR IS 1 ROW BACK
      B=DZERO
      BHNEW=DZERO
      IF(I.NE.1) THEN
        B=CC(NRB)
        BHNEW=B*(HNEW(NRL)-HNEW(N))
      ENDIF
C
C-------NEIGHBOR IS 1 ROW AHEAD
      H=DZERO
      HHNEW=DZERO
      IF(I.NE.NROW) THEN
        H=CC(NRH)
        HHNEW=H*(HNEW(NRN)-HNEW(N))
      ENDIF
C
C-------NEIGHBOR IS 1 COLUMN BACK
      D=DZERO
      DHNEW=DZERO
      IF(J.NE.1) THEN
        D=CR(NCD)
        DHNEW=D*(HNEW(NCL)-HNEW(N))
      ENDIF
C
C-------NEIGHBOR IS 1 COLUMN AHEAD
      F=DZERO
      FHNEW=DZERO
      IF(J.NE.NCOL) THEN
        F=CR(NCF)
        FHNEW=F*(HNEW(NCN)-HNEW(N))
      ENDIF
C
C-------NEIGHBOR IS 1 LAYER BEHIND
      Z=DZERO
      ZHNEW=DZERO
      IF(K.NE.1) THEN
        Z=CV(NLZ)
        ZHNEW=Z*(HNEW(NLL)-HNEW(N))
      ENDIF
C
C-------NEIGHBOR IS 1 LAYER AHEAD
      S=DZERO
      SHNEW=DZERO
      IF(K.NE.NLAY) THEN
        S=CV(NLS)
        SHNEW=S*(HNEW(NLN)-HNEW(N))
      ENDIF
C
      IF(I.EQ.NROW) CC(N)=0.
      IF(J.EQ.NCOL) CR(N)=0.
C-------15JUN1993 SKIP CALCULATIONS AND MAKE CELL INACTIVE IF ALL
C                 SURROUNDING CELLS ARE INACTIVE
      IF(B+H+D+F+Z+S.EQ.0.) THEN
        IBOUND(N)=0
        HCOF(N)=0.
        RHS(N)=0.
        GO TO 115
      ENDIF
C
C-------CALCULATE THE RESIDUAL AND STORE IT IN RHS.  TO SCALE A,
C-------CALCULATE THE DIAGONAL OF THE A MATRIX, AND STORE IT IN HCOF.
      E=-Z-B-D-F-H-S
      RRHS=RHS(N)
      HHCOF=HNEW(N)*HCOF(N)
      RHS(N)=RRHS-ZHNEW-BHNEW-DHNEW-HHCOF-FHNEW-HHNEW-SHNEW
      IF(NORM.EQ.1) HCOF(N)=HCOF(N)+E
      IF(IBOUND(N).LT.0) RHS(N)=0.
C-------ADDED FOR SENSITIVITY CALCULATIONS 9/1/91
      IF(IU.NE.0.AND.IP.GT.0) THEN
        IF(I.EQ.1.AND.J.EQ.1.AND.K.EQ.1) THEN
          SN=0.
          SP=0.
          SR=0.
        ENDIF
        SR=SR+RHS(N)
        IF(RRHS.LT.0.) SN=SN+RRHS
        IF(RRHS.GT.0.) SP=SP+RRHS
        IF(-ZHNEW.LT.0.) SN=SN-ZHNEW
        IF(-ZHNEW.GT.0.) SP=SP-ZHNEW
        IF(-BHNEW.LT.0.) SN=SN-BHNEW
        IF(-BHNEW.GT.0.) SP=SP-BHNEW
        IF(-DHNEW.LT.0.) SN=SN-DHNEW
        IF(-DHNEW.GT.0.) SP=SP-DHNEW
        IF(-HHCOF.LT.0.) SN=SN-HHCOF
        IF(-HHCOF.GT.0.) SP=SP-HHCOF
        IF(-FHNEW.LT.0.) SN=SN-FHNEW
        IF(-FHNEW.GT.0.) SP=SP-FHNEW
        IF(-HHNEW.LT.0.) SN=SN-HHNEW
        IF(-HHNEW.GT.0.) SP=SP-HHNEW
        IF(-SHNEW.LT.0.) SN=SN-SHNEW
        IF(-SHNEW.GT.0.) SP=SP-SHNEW
      ENDIF
  115 CONTINUE
C
C-------SCALE CC,CR,CV,RHS AND HNEW IF NORM=1.
      IF(NORM.EQ.1) THEN
        DO 120 K=1,NLAY
        DO 120 I=1,NROW
        DO 120 J=1,NCOL
          N=J+(I-1)*NCOL+(K-1)*NRC
          IF(IBOUND(N).EQ.0) GO TO 120
          HHCOF=SQRT(-HCOF(N))
          IF(N.LE.(NODES-NCOL).AND.CC(N).GT.0.)
     #      CC(N)=CC(N)/(HHCOF*(SQRT(-HCOF(N+NCOL))))
          IF(CR(N).GT.0.) CR(N)=CR(N)/(HHCOF*(SQRT(-HCOF(N+1))))
          IF(N.LE.(NODES-NRC).AND.CV(N).GT.0.)
     #      CV(N)=CV(N)/(HHCOF*(SQRT(-HCOF(N+NRC))))
          HNEW(N)=HNEW(N)*HHCOF
          RHS(N)=RHS(N)/HHCOF
  120   CONTINUE
      ENDIF
C
C-------CALCULATE PARAMETER B OF THE POLYNOMIAL PRECONDITIONING METHOD
      IF(NPCOND.NE.2) GO TO 152
      IF(NBPOL.EQ.2) THEN
        BPOLY=2
        GO TO 151
      ENDIF
      DO 150 K=1,NLAY
      DO 150 I=1,NROW
      DO 150 J=1,NCOL
C
      N=J+(I-1)*NCOL+(K-1)*NRC
      IF(IBOUND(N).LE.0)GO TO 150
C
      NCF=N
      NCD=N-1
      NRB=N-NCOL
      NRH=N
      NLS=N
      NLZ=N-NRC
C
      B=DZERO
      IF(I.NE.1) B=CC(NRB)
      H=DZERO
      IF(I.NE.NROW) H=CC(NRH)
      D=DZERO
      IF(J.NE.1) D=CR(NCD)
      F=DZERO
      IF(J.NE.NCOL) F=CR(NCF)
      Z=DZERO
      IF(K.NE.1) Z=CV(NLZ)
      S=DZERO
      IF(K.NE.NLAY) S=CV(NLS)
C
C-------NOTE : ABS. VAL. OF THE DIAG. OF THE SCALED A MATRIX IS 1.
      HHCOF=HCOF(N)
      IF(NORM.EQ.1) HHCOF=DONE
      T=DABS(Z)+DABS(B)+DABS(D)+ABS(HHCOF)+DABS(F)+DABS(H)+DABS(S)
      IF(T.GT.BPOLY) BPOLY=T
  150 CONTINUE
  151 CONTINUE
C
C-------CALCULATE ITERATION PARAMETERS FOR POLYNOMIAL PRECONDITIONING
C-------METHOD FOR A NEGATIVE DEFINITE MATRIX.
        C0=(15./32.)*(BPOLY**3)
        C1=(27./16.)*(BPOLY**2)
        C2=(9./4.)*BPOLY
  152 CONTINUE
C
C-------START INTERNAL ITERATIONS
      IITER=0
      IF(KITER.EQ.1) NITER=0
      ICNVG=0
      IICNVG=0
153   CONTINUE
      IITER=IITER+1
      NITER=NITER+1
C
C-------INITIALIZE VARIABLES THAT TRACK MAXIMUM HEAD CHANGE AND RESIDUAL
C-------VALUE DURING EACH ITERATIONS
      BIGH=0.
      BIGR=0.
C-------INITIALIZE DEL (ADDED 01MAY1993)
      DEL=0.
C
C
C-------CHECK NPCOND FOR PRECONDITIONING TYPE AND EXECUTE PROPER CODE
      IF(NPCOND.EQ.2) GO TO 165
C
C-------CHOLESKY PRECONDITIONING
C
C-------STEP THROUGH CELLS TO CALCULATE THE DIAGONAL OF THE CHOLESKY
C-------MATRIX (FIRST INTERNAL ITERATION ONLY) AND THE INTERMEDIATE
C-------SOLUTION.  STORE THEM IN CD AND V, RESPECTIVELY.
154   DO 155 K=1,NLAY
      DO 155 I=1,NROW
      DO 155 J=1,NCOL
C
      N=J+(I-1)*NCOL+(K-1)*NRC
      IF(IBOUND(N).LE.0)GO TO 155
C
C-------CALCULATE V
      H=DZERO
      VCC=DZERO
      IC=N-NCOL
      IF(I.NE.1) THEN
        H=CC(IC)
        IF(CD(IC).NE.0.) VCC=H*V(IC)/CD(IC)
      ENDIF
C
      F=DZERO
      VCR=DZERO
      IR=N-1
      IF(J.NE.1) THEN
        F=CR(IR)
        IF(CD(IR).NE.0.) VCR=F*V(IR)/CD(IR)
      ENDIF
C
      S=DZERO
      VCV=DZERO
      IL=N-NRC
      IF(K.NE.1) THEN
        S=CV(IL)
        IF(CD(IL).NE.0.) VCV=S*V(IL)/CD(IL)
      ENDIF
      V(N)=RHS(N)-VCR-VCC-VCV
C
C-------CALCULATE CD - FIRST INTERNAL ITERATION ONLY
      IF(IITER.EQ.1) THEN
          CDCR=DZERO
          CDCC=DZERO
          CDCV=DZERO
          FCC=DZERO
          FCR=DZERO
          FCV=DZERO
          IF(IR.GT.0.AND.CD(IR).NE.0.) CDCR=(F**2)/CD(IR)
          IF(IC.GT.0.AND.CD(IC).NE.0.) CDCC=(H**2)/CD(IC)
          IF(IL.GT.0.AND.CD(IL).NE.0.) CDCV=(S**2)/CD(IL)
          IF(NPCOND.EQ.1) THEN
              IF(IR.GT.0) THEN
                  FV=CV(IR)
C                 MODIFIED FROM HILL(1990) 9/27/90: 2 REPLACES 1
                  IF(K.EQ.NLAY.AND.((J+I).GT.2)) FV=DZERO
                  IF(CD(IR).NE.0.) FCR=(F/CD(IR))*(CC(IR)+FV)
              ENDIF
              IF(IC.GT.0) THEN
                  FV=CV(IC)
                  IF(K.EQ.NLAY.AND.(I.GT.1)) FV=DZERO
                  IF(CD(IC).NE.0.) FCC=(H/CD(IC))*(CR(IC)+FV)
              ENDIF
              IF(IL.GT.0) THEN
                  IF(CD(IL).NE.0.) FCV=(S/CD(IL))*(CR(IL)+CC(IL))
              ENDIF
          ENDIF
          IF(NORM.EQ.0) THEN
             B=DZERO
             H=DZERO
             D=DZERO
             F=DZERO
             Z=DZERO
             S=DZERO
             IF(I.NE.1) B=CC(IC)
             IF(I.NE.NROW) H=CC(N)
             IF(J.NE.1) D=CR(IR)
             IF(J.NE.NCOL) F=CR(N)
             IF(K.NE.1) Z=CV(IL)
             IF(K.NE.NLAY) S=CV(N)
             HHCOF=HCOF(N)-Z-B-D-F-H-S
          ENDIF
          IF(NORM.EQ.1) HHCOF=-DONE
          CD(N)=(DONE+DEL)*HHCOF-CDCR-CDCC-CDCV-RELAX*(FCR+FCC+FCV)
         IF(CD1.EQ.0..AND.CD(N).NE.0.) CD1=CD(N)
C--------.LT. CHANGED TO .LE. 01SEPT1991
         IF(CD(N)*CD1.LE.0.) THEN
C--------CHANGED 510 FORMAT 01SEPT1991 AND 20MAR1992
C--------CHANGED 510 FORMAT AND WRITE STATEMENT AND ADDED DEL 01MAY1993
           DEL=1.5*DEL+.001
           IF(DEL.GT..5) THEN
             WRITE(IOUT,510)
             STOP
           ENDIF
           GO TO 154
         ENDIF
      ENDIF
510   FORMAT(' MATRIX IS SEVERELY NON-DAGONALLY DOMINANT.  CHECK INPUT',
     1     ' FILES.  STOP EXECUTION')
C
  155 CONTINUE
C
C-------STEP THROUGH EACH CELL AND SOLVE FOR S OF THE CONJUGATE
C-------GRADIENT ALGORITHM BY BACK SUBSTITUTION. STORE RESULT IN SS.
      DO 160 KK=NLAY,1,-1
      DO 160 II=NROW,1,-1
      DO 160 JJ=NCOL,1,-1
C
      N=JJ+(II-1)*NCOL+(KK-1)*NRC
      IF(IBOUND(N).LE.0)GO TO 160
C
      NC=N+1
      NR=N+NCOL
      NL=N+NRC
      SSCR=DZERO
      SSCC=DZERO
      SSCV=DZERO
      IF(JJ.NE.NCOL) SSCR=CR(N)*SS(NC)/CD(N)
      IF(II.NE.NROW) SSCC=CC(N)*SS(NR)/CD(N)
      IF(KK.NE.NLAY) SSCV=CV(N)*SS(NL)/CD(N)
      VN=V(N)/CD(N)
      SS(N)=VN-SSCR-SSCC-SSCV
  160 CONTINUE
C-------SKIP OVER OTHER PRECONDITIONING TYPES
      GO TO 199
  165 CONTINUE
C
C-------POLYNOMIAL PRECONDITIONING
      DO 170 N=1,NODES
      V(N)=RHS(N)
  170 CONTINUE
      CALL SPCG2E(IBOUND,RHS,HCOF,CR,CC,CV,V,SS,C2,NORM,NCOL,NROW,
     1            NLAY,NODES)
      CALL SPCG2E(IBOUND,RHS,HCOF,CR,CC,CV,SS,V,C1,NORM,NCOL,NROW,
     1            NLAY,NODES)
      CALL SPCG2E(IBOUND,RHS,HCOF,CR,CC,CV,V,SS,C0,NORM,NCOL,NROW,
     1            NLAY,NODES)
  199 CONTINUE
C
C-------CALCULATE P OF THE CONJUGATE GRADIENT ALGORITHM
      SROLD=SRNEW
      SRNEW=DZERO
      DO 200 N=1,NODES
        IF(IBOUND(N).LE.0)GO TO 200
        SRNEW=SRNEW+SS(N)*RHS(N)
  200 CONTINUE
C
      IF(IITER.EQ.1) THEN
         DO 205 N=1,NODES
  205    P(N)=SS(N)
      ELSE
         DO 210 N=1,NODES
  210    P(N)=SS(N)+(SRNEW/SROLD)*P(N)
      ENDIF
C
C-------CALCULATE ALPHA OF THE CONJUGATE GRADIENT ROUTINE.
C-------FOR THE DENOMINATOR OF ALPHA, MULTIPLY THE MATRIX A BY THE
C-------VECTOR P, AND STORE IN V; THEN MULTIPLY P BY V.  STORE IN PAP.
      PAP=DZERO
      DO 290 K=1,NLAY
      DO 290 I=1,NROW
      DO 290 J=1,NCOL
C
      N=J+(I-1)*NCOL+(K-1)*NRC
      V(N)=0.
      IF(IBOUND(N).LE.0)GO TO 290
C
      NRN=N+NCOL
      NRL=N-NCOL
      NCN=N+1
      NCL=N-1
      NLN=N+NRC
      NLL=N-NRC
C
      NCF=N
      NCD=NCL
      NRB=NRL
      NRH=N
      NLS=N
      NLZ=NLL
C
      B=DZERO
      IF(I.NE.1) B=CC(NRB)
      H=DZERO
      IF(I.NE.NROW)H=CC(NRH)
      D=DZERO
      IF(J.NE.1) D=CR(NCD)
      F=DZERO
      IF(J.NE.NCOL) F=CR(NCF)
      Z=DZERO
      IF(K.NE.1) Z=CV(NLZ)
      S=DZERO
      IF(K.NE.NLAY) S=CV(NLS)
C
      IF(NORM.EQ.0) PN=P(N)
      IF(NORM.EQ.1) PN=DZERO
      BHNEW=DZERO
      HHNEW=DZERO
      DHNEW=DZERO
      FHNEW=DZERO
      ZHNEW=DZERO
      SHNEW=DZERO
      IF(NRL.GT.0) BHNEW=B*(P(NRL)-PN)
      IF(NRN.LE.NODES) HHNEW=H*(P(NRN)-PN)
      IF(NCL.GT.0) DHNEW=D*(P(NCL)-PN)
      IF(NCN.LE.NODES) FHNEW=F*(P(NCN)-PN)
      IF(NLL.GT.0) ZHNEW=Z*(P(NLL)-PN)
      IF(NLN.LE.NODES) SHNEW=S*(P(NLN)-PN)
C
C-------CALCULATE THE PRODUCT OF MATRIX A AND VECTOR P AND STORE
C-------RESULT IN V.
      PN=HCOF(N)*P(N)
      IF(NORM.EQ.1) PN=-P(N)
      VN=ZHNEW+BHNEW+DHNEW+PN+FHNEW+HHNEW+SHNEW
      V(N)=VN
      PAP=PAP+P(N)*VN
  290 CONTINUE
C
C-------CALCULATE ALPHA
      ALPHA=1.
      IF(PAP.EQ.0..AND.MXITER.EQ.1) THEN
        WRITE(IOUT,520)
        STOP
      ENDIF
520     FORMAT(/,' CONJUGATE-GRADIENT METHOD FAILED.',/,' SET MXITER',
     1           ' GREATER THAN ONE AND TRY AGAIN.  STOP EXECUTION')
      IF(PAP.NE.0.) ALPHA=SRNEW/PAP
C
C-------CALCULATE NEW HEADS AND RESIDUALS, AND SAVE THE LARGEST
C-------CHANGE IN HEAD AND THE LARGEST VALUE OF THE RESIDUAL.
      DO 300 K=1,NLAY
      DO 300 I=1,NROW
      DO 300 J=1,NCOL
C
      N=J+(I-1)*NCOL+(K-1)*NRC
      IF(IBOUND(N).LE.0) GO TO 300
C
C-------HEAD
      HCHGN=ALPHA*P(N)
      IF(DABS(HCHGN).GT.ABS(BIGH)) THEN
          BIGH=HCHGN
          IH=I
          JH=J
          KH=K
          NH=N
      ENDIF
      HNEW(N)=HNEW(N)+DAMP*HCHGN
C
C--------RESIDUAL (V IS THE PRODUCT OF MATRIX A AND VECTOR P)
      RCHGN=-ALPHA*V(N)
      RHS(N)=RHS(N)+DAMP*RCHGN
      IF(ABS(RHS(N)).GT.ABS(BIGR)) THEN
          BIGR=RHS(N)
          IR=I
          JR=J
          KR=K
          NR=N
      ENDIF
  300 CONTINUE
C
C-------UNSCALE LARGEST CHANGE IN HEAD AND LARGEST RESIDUAL, AND
C-------CHECK THE CONVERGENCE CRITERION
       IF(NORM.EQ.1) THEN
          BIGH=BIGH/SQRT(-HCOF(NH))
          BIGR=BIGR*SQRT(-HCOF(NR))
       ENDIF
       IF(MXITER.EQ.1) THEN
          IF(ABS(BIGH).LE.HCLOSE.AND.ABS(BIGR).LE.RCLOSE) ICNVG=1
       ELSE
          IF(IITER.EQ.1.AND.
     1       ABS(BIGH).LE.HCLOSE.AND.ABS(BIGR).LE.RCLOSE) ICNVG=1
       ENDIF
       IF(ABS(BIGH).LE.HCLOSE.AND.ABS(BIGR).LE.RCLOSE) IICNVG=1
C
C-------STORE THE LARGEST UNSCALED HEAD CHANGE AND RESIDUAL VALUE
C-------(THIS ITERATION) AND THEIR LOCATIONS.
      II=NITER
      HCHG(II)=BIGH
      LHCH(1,II)=KH
      LHCH(2,II)=IH
      LHCH(3,II)=JH
C
      RCHG(II)=BIGR
      LRCH(1,II)=KR
      LRCH(2,II)=IR
      LRCH(3,II)=JR
C
      IT1(II)=0
      IF(IITER.EQ.1) IT1(II)=1
C-------GO TO NEXT INTERNAL ITERATION IF CONVERGENCE HAS NOT BEEN
C-------REACHED AND IITER IS LESS THAN ITER1
      IF(MXITER.EQ.1) THEN
        IF(ICNVG.EQ.0.AND.IITER.LT.ITER1) GO TO 153
      ELSE
        IF(IICNVG.EQ.0.AND.IITER.LT.ITER1) GO TO 153
      ENDIF
C
C-------UNSCALE CR,CC,CV AND HNEW
      IF(NORM.EQ.1) THEN
      DO 310 N=1,NODES
      IF(IBOUND(N).EQ.0) GO TO 310
      HHCOF=SQRT(-HCOF(N))
      IF(N.LE.(NODES-NCOL).AND.CC(N).GT.0.)
     #   CC(N)=CC(N)*(HHCOF*(SQRT(-HCOF(N+NCOL))))
      IF(N.LE.(NODES-1).AND.CR(N).GT.0.)
     #   CR(N)=CR(N)*(HHCOF*(SQRT(-HCOF(N+1))))
      IF(N.LE.(NODES-NRC).AND.CV(N).GT.0.)
     #   CV(N)=CV(N)*(HHCOF*(SQRT(-HCOF(N+NRC))))
      HNEW(N)=HNEW(N)/HHCOF
  310 CONTINUE
      ENDIF
C
C-------IF END OF TIME STEP, PRINT # OF ITERATIONS THIS STEP
      IF(ICNVG.EQ.0 .AND. KITER.NE.MXITER) GO TO 600
      IF(MUTPCG.GT.1) GO TO 590
      IF(KSTP.EQ.1) WRITE(IOUT,500)
  500 FORMAT(1H0)
      WRITE(IOUT,501) KITER,KSTP,KPER,NITER
  501 FORMAT(1X,I5,' CALLS TO PCG ROUTINE FOR TIME STEP',I4,
     1' IN STRESS PERIOD',I3,/1X,I5,' TOTAL ITERATIONS')
      IF(MUTPCG.EQ.1) GO TO 590
C
C-------PRINT HEAD CHANGE EACH ITERATION IF PRINTOUT INTERVAL IS REACHED
      IF(ICNVG.EQ.0 .OR. KSTP.EQ.NSTP .OR. MOD(KSTP,IPRPCG).EQ.0)
     1      CALL SPCG2P(HCHG,LHCH,RCHG,LRCH,
     2                  ITER1,NITER,MXITER,IOUT,NPCOND,BPOLY,IT1)
590   NITER=0
C
C-------RETURN
600   RETURN
C
      END
      SUBROUTINE SPCG2P(HCHG,LHCH,RCHG,LRCH,ITER1,
     1                  NITER,MXITER,IOUT,NPCOND,BPOLY,IT1)
C
C
C-----VERSION 0001 01MAY1988 SPCG2P
C     ******************************************************************
C     PRINT MAXIMUM HEAD CHANGE AND RESIDUAL VALUE FOR EACH ITERATION
C                           DURING A TIME STEP
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
C
      DIMENSION HCHG(MXITER*ITER1), LHCH(3,MXITER*ITER1)
      DIMENSION RCHG(MXITER*ITER1), LRCH(3,MXITER*ITER1)
      DIMENSION IT1(MXITER*ITER1)
C     ------------------------------------------------------------------
C
      IF(NPCOND.EQ.2) WRITE(IOUT,2) BPOLY
    2 FORMAT(1H0,'B OF THE POLYNOMIAL PRECONDITIONING METHOD: ',E12.4)
      WRITE(IOUT,5)
    5 FORMAT(1H0,'MAXIMUM HEAD CHANGE FOR EACH ITERATION (1 INDICATES ',
     A      'THE FIRST INNER ITERATION):'/
     1    1H0,4('    HEAD CHANGE  LAYER,ROW,COL')/1X,120('-'))
      WRITE (IOUT,10) (IT1(J),HCHG(J),(LHCH(I,J),I=1,3),J=1,NITER)
      WRITE(IOUT,11)
      WRITE(IOUT,15)
   15 FORMAT(1H0,'MAXIMUM RESIDUAL FOR EACH ITERATION (1 INDICATES ',
     A      'THE FIRST INNER ITERATION):'/
     1    1H0,4('     RESIDUAL    LAYER,ROW,COL')/1X,120('-'))
      WRITE (IOUT,10) (IT1(J),RCHG(J),(LRCH(I,J),I=1,3),J=1,NITER)
      WRITE(IOUT,11)
C
      RETURN
C
   10 FORMAT((1X,4(2X,I1,1X,G12.4,' (',I3,',',I3,',',I3,')')))
   11 FORMAT(1H0)
C
      END
      SUBROUTINE SPCG2E(IBOUND,RHS,HCOF,CR,CC,CV,VIN,VOUT,C,NORM,NCOL,
     1                  NROW,NLAY,NODES)
C
C
C-----VERSION 0001 01MAY1989 SPCG2E
C     ******************************************************************
C           MATRIX MULTIPLICATIONS FOR POLYNOMIAL PRECONDITIONING
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
C
      DOUBLE PRECISION VN,CRHS,Z,B,D,F,H,S,ZV,BV,DV,FV,HV,SV,DZERO
      DIMENSION IBOUND(NODES),CR(NODES),CC(NODES),CV(NODES),
     1  RHS(NODES),VIN(NODES),VOUT(NODES),HCOF(NODES)
C     ------------------------------------------------------------------
C
      DZERO=0.
      NRC=NROW*NCOL
      DO 290 K=1,NLAY
      DO 290 I=1,NROW
      DO 290 J=1,NCOL
C
      N=J+(I-1)*NCOL+(K-1)*NRC
      VOUT(N)=0.
      IF(IBOUND(N).LE.0)GO TO 290
C
      NRN=N+NCOL
      NRL=N-NCOL
      NCN=N+1
      NCL=N-1
      NLN=N+NRC
      NLL=N-NRC
C
      NCF=N
      NCD=NCL
      NRB=NRL
      NRH=N
      NLS=N
      NLZ=NLL
C
      B=DZERO
      BV=DZERO
      IF(I.NE.1.AND.IBOUND(NRL).GE.0) THEN
        B=CC(NRB)
        BV=B*VIN(NRL)
      ENDIF
      H=DZERO
      HV=DZERO
      IF(I.NE.NROW.AND.IBOUND(NRN).GE.0) THEN
        H=CC(NRH)
        HV=H*VIN(NRN)
      ENDIF
      D=DZERO
      DV=DZERO
      IF(J.NE.1.AND.IBOUND(NCL).GE.0) THEN
        D=CR(NCD)
        DV=D*VIN(NCL)
      ENDIF
      F=DZERO
      FV=DZERO
      IF(J.NE.NCOL.AND.IBOUND(NCN).GE.0) THEN
        F=CR(NCF)
        FV=F*VIN(NCN)
      ENDIF
      Z=DZERO
      ZV=DZERO
C      IF STATEMENT REARRANGED 01JUN1993
      IF(K.NE.1) THEN
        IF(IBOUND(NLL).GE.0) THEN
          Z=CV(NLZ)
          ZV=Z*VIN(NLL)
        ENDIF
      ENDIF
      S=DZERO
      SV=DZERO
      IF(K.NE.NLAY.AND.IBOUND(NLN).GE.0) THEN
        S=CV(NLS)
        SV=S*VIN(NLN)
      ENDIF
C
C-------CALCULATE THE PRODUCT OF MATRIX A AND VECTOR VIN AND STORE
C------ RESULT IN VOUT
      VN=HCOF(N)*VIN(N)
      IF(NORM.EQ.1) VN=-VIN(N)
      CRHS=C*RHS(N)
      VOUT(N)=CRHS+ZV+BV+DV+VN+FV+HV+SV
  290 CONTINUE
      RETURN
      END
